Self-organized and directed branching results in optimal coverage in developing dermal lymphatic networks

Branching morphogenesis is a ubiquitous process that gives rise to high exchange surfaces in the vasculature and epithelial organs. Lymphatic capillaries form branched networks, which play a key role in the circulation of tissue fluid and immune cells. Although mouse models and correlative patient data indicate that the lymphatic capillary density directly correlates with functional output, i.e., tissue fluid drainage and trafficking efficiency of dendritic cells, the mechanisms ensuring efficient tissue coverage remain poorly understood. Here, we use the mouse ear pinna lymphatic vessel network as a model system and combine lineage-tracing, genetic perturbations, whole-organ reconstructions and theoretical modeling to show that the dermal lymphatic capillaries tile space in an optimal, space-filling manner. This coverage is achieved by two complementary mechanisms: initial tissue invasion provides a non-optimal global scaffold via self-organized branching morphogenesis, while VEGF-C dependent side-branching from existing capillaries rapidly optimizes local coverage by directionally targeting low-density regions. With these two ingredients, we show that a minimal biophysical model can reproduce quantitatively whole-network reconstructions, across development and perturbations. Our results show that lymphatic capillary networks can exploit local self-organizing mechanisms to achieve tissue-scale optimization.

as a magnified image (n=6 ear pinna, representing 6 mice).Branch points are indicated with orange nodes and tips with green nodes.Quantification of the network parameters are shown in Fig. 1D-F Welch's t-test was used for measuring statistical significance.p=0.0003 (P8 vs. P16), p=0.0004 (P8 vs. P21), and p=0.93 (P16 vs. P21).I) Fluctuation exponents of the lymphatic capillary (LYVE1+) or total LV network (VEGFR3) calculated from P21 ventral networks stained with anti-LYVE1 and anti-VEGFR3 (n=10 ear pinna, representing 9 mice each), showing that both stainings display minimal spatial fluctuations.p=0.49, two-sided Welch's t-test.Data are shown as mean values +/-SD.J) Examples of density fluctuation spectra for a P13 and a P21 ventral networks (dots) and power-law fit (taken from x=10 to 750, lines).Scale bars in A) 1000µm,D-F) 500µm, and G) 1000µm or 500µm (zoom-in).Source data for Supplementary Fig. 1H-J are provided as a Source Data file.
A) P6 and P21 hair follicles (green) in ventral mouse ear pinna as highlighted by Sox9 promoter-driven EGFP.Scale bars are 200µm.B-C) Dot blot shows the mean +/-SD of median distances of hair follicles in each analyzed ear pinna at P6 and P21 (p<0.0001),whereas the scatter blots show the distances of all the hair follicles to their nearest neighbour (x-axes) as a function of distance to the outer border of the ear pinna (y-axes).n=7 ear pinna representing 5 mice (P6) or n=4 ear pinna representing 3 mice (P21).In B) two-sided Welch's t-test was used for measuring statistical significance.D) The flattened overview image shows anti-LYVE1 stained (red) superficial LVs invading the ventral dermis at the base of the P7 ear pinna.The zoom-in images show a single optical section of LYVE1+ LVs and EdU (white) that incorporated the DNA upon 4h pulse labeling.The zoom-in image framed in white highlights a tip of an LV tree whereas the zoom-in image, below, framed in yellow highlights a part of the trunk of an branching (all branches born from side-branching events in red, all others in black), we find that the system becomes rather insensitive to the value of the tip-branching rate rb, consistent with side-branching being able to rescue low branching rates (G, left to right) -see also Fig. 2E,G.I-K) Effect of different values of side-branching rate (different colors) on the final network morphology (giving rise to larger densities), as well as on the temporal evolution of different morphometric parameters.In particular, branch length is expected to linearly increase through time in the absence of side-branching due to overall ear dilation, while the effect can be dampened or even reversed as a function of the side-branching rate giving rise to smaller branches (I).As expected, increasing values of side-branching rate in the second phase of branching also result in increasing branch numbers (J) and total network length (K).

L,M)
Value of fluctuation exponents in simulations under different theoretical assumptions (corresponding to respectively Fig. 2E,G).Fluctuation exponents do not converge towards 0.5 (minimal spatial fluctuations) upon an increase of tip-branching rate (L, calculated from final P21 snapshots), but do rapidly converge towards 0.5 as side-branching proceeds (M, calculated as a function of time/number of side-branching events).Source data for Supplementary Fig. 3L-M  Sensitivity of the model prediction to changes in the initial conditions, where we have considered a sixth initial tree growing from the ear tip (left).This produces networks that are highly similar to those shown in Fig. 2, as shown quantitatively by similar simulated clone size distribution (right panel, compare with Fig. 3F, see also Supplementary Movie 12).Source data for Supplementary Fig. 6A-B shows the mean +/-SD percentage of side-branching events (sprouting of an existing vessel) of all branching events (side-branching + tip-splitting) at P6 (n=5 ear pinna and mice), P8 (n=5 ear pinna and mice), and P13 (n=5 ear pinna, representing 3 mice).Two-sided Welch's t-test was used for measuring statistical significance.p=0.04 (P6 vs. P8), p=0.001 (P6 vs. P13), p=0.04 (P8 vs. P13).These samples were part of the data set also shown in Fig. 1B, Supplementary Fig. 1C-F and H and Supplementary Movies 2-6.C) Schematic illustration of the vertex configuration factor σ which quantifies the angular configurations around branching events by reweighting the ratio of the largest angle α to π with the ratio of the two smaller angles γ to β.The drawings illustrate values of σ for different configurations: σ=1 for perfectly lateral branching modes (left), σ=2/3 for symmetric configurations (middle), and σ<2/3 for "forklike" bifurcations (right).D) Fraction of branching configurations with σ>2/3 relative to all calculated σ values from entire LV networks across different developmental stages.A larger fraction indicates an increasing occurrence of symmetric and lateral branching modes as illustrated in panel C.The column "P13 (side)" shows vertex configuration factors σ for manually identified side-branching events (n=236) from P13 datasets.E) Cumulative histograms of calculated vertex configuration factors σ across different developmental stages, and for the manually identified P13 side-branches (rightmost panel).Dotted (orange) and dashed (red) vertical lines indicate the mean of samples and the value σ=2/3, respectively.
Source data for Supplementary Fig. 7B and D-E  Supplementary Theory -Self-organized and directed branching results in optimal coverage in developing dermal lymphatic networks In this Supplemental Theory Note, we provide additional details on the modelling approach, as well as sensitivity analyses on how di↵erent models and parameters a↵ect the structure of in silico lymphatic networks compared to our in vivo dataset.
1 Models of lymphatic network formation as a branching and annihilating random walk Our data shown in Figure 1 supports that post-natal lymphatic network formation in the mouse ear pinna dermis proceeds via two complementary processes: a first phase of invasive branching morphogenesis (from P4 to P10-11), where a small number of tree-like networks invade the ear via iterative rounds of branching and elongation, as well as a second phase (P11-P21), where the ear is already filled by lymphatic networks, but where we still observe an increase in the total number of branches in a given ear (Figure 1).

Description of the model of branching and annihilating random walks
Given that the first phase of branching morphogenesis gives rise to networks which are highly stochastic from one mouse to the next, as well as ine cient at space-filling (Figure 2), we first consider the simplest self-organized model for this process as a branching and annihilating random walk.This type of models, which has been proposed to understand branched organ morphogenesis (such as mammary gland, kidney or pancreas) as well as neuronal morphology [1], considers that branch structures emerge via a set of simple rules.Actively growing tips are responsible for the bulk of the organ growth, and drive local ductal elongation by undergoing a persistent random walk (note that we do not model the details of the mechanics involved in this, which presumably requires both cell migration and proliferation).Furthermore, a given tip has at any time point a given probability to branch into two tips (branching random walk).Finally, these two rules alone would give rise to an uncontrolled exponential growth of the structure, so that we need to posit a mechanism for tip termination/annihilation.Interestingly, ear lymphatic networks form largely in 2D, with very little overlap between vessels, and in particular show consistently tips that are arrested in growth at a given, small distance, from another vessel.This is close to previous observation in organs such as mammary gland, and can be modelled by density-dependent tip termination, where a given active tip irreversibly becomes inactive (i.e.cannot elongate or branch anymore) when it gets too close to a duct.These three rules have been shown to be enough to give rise to self-organized branched organ growth at a set density [2].
More precisely, considering active tips as particles A and inactive vessels left behind by tip elongation or terminal as particles I, the model has the following parameters (note that for simplicity, we consider all tips to have the same underlying parameters): • tips elongate at an instantaneous speed v 0 , with angular di↵usion D r (persistent random walk), giving rise to vessels in their wake.(A !A + I process) • tips branch into two tips probabilistically at rate r b (A !A + A process).Note that this means that the branch length will be exponentially distributed, with characteristic length scale related to v 0 /r b .
• tips terminate irreversibly if they are within a radius R a of other vessels.(A !0 process) We perform stochastic numerical simulations of this process, where at every time interval dt, a given tip i moves by a distance v 0 dt in its preferred direction q i .This angle q i changes gradually at every time step (persistent random walk), by the addition of a small random component picked from a uniform distribution between [ d q , d q ], with d q = p/10.At every time point, we pick a random number and compare it to r b dt to determine whether this tip also branches into two tips.If so, we terminate the tip i and create two new tips with angular direction q i + q 1 and q i q 2 , where q 1 , q 2 are random branching angles picked from a uniform distribution [p/6, 2p/3].Finally, we check at every time point whether the tip is below a distance R a of any other vessels, or outside the boundary of the simulations, in which cases we terminate its growth.
Note that in two-dimensions, the exact value of R a has a small influence on the network growth (as long as it is small compared to other length scales) as the probability of two branches to intersect is 1.The exact value of the rotational di↵usivity and branch angle distribution also has a mild e↵ect on the overall dynamics [2].
An important feature of this model is that termination/annihilation is irreversible: once a given branched region has no active tips anymore, it is "frozen".This manifests into very large (non-equilibrium) spatial fluctuations, as quantified in analogy to active matter systems from statistical physics, by quantifying the spatial standard deviation vs spatial averages of network density at multiple length scales L.More specifically, we divide a given system into boxes of size LxL, and count for each box the network density r L .We then calculate population-level standard deviation and averages across all boxes (s (r) L and < r > L ), and plot these two quantifies as a function of one another across many values of L. Equilibrium physics imposes that s (r) L µ (< r > L ) a , with a = 0.5, whereas exponents a > 0.5 indicate giant fluctuations (note that a must be small than 1, as this corresponds to maximal standard deviation, i.e. all the network to be in one box).
Interestingly, calculating the fluctuation exponent a for the P13 data revealed significant departure from minimal fluctuations (Figure 2A,B), with values of a = 0.6 which were fully consistent with previous exponents found for our model of branching and annihilating random walk, across di↵erent parameter regimes (Figure 2D-E, Supplementary Figure 3A-F).This qualitatively suggests that this model might be able to capture some features of our data.

Incorporation of tip-vessel repulsion in the model
In addition to the simple rules for branching and annihilating random walk discussed above, a possible extension of the model, that has been studied in particular in the context of branched neurons, is that tips do not simply terminate in the vicinity of surrounding branches, but also actively avoid them, by processes such as adhesion-mediated repulsion [3] or local gradient sensing [4].In our framework, this can be captured by considering that at every time point, a given tip locally senses vessels density in a given radius R rep and adapts its direction of growth accordingly.
More specifically, for a tip i with direction of growth q i and position ri , we compute the local density vector R = Â j (r i r j ) Rrep , where index j lists all the neighbors of i (for computational convenience, we only calculate it for |r i r j | < 2R rep given that vessels further away have vanishing contributions).This local density vector R is defined such that it points towards regions of low network density (with density being "measured" locally as vessels further contributing less compared to vessels nearby the tip).Then, in the elongation step, instead of updating as before tip coordinates (x i , y i ) as (x i + v 0 dtcos(q i ), y i + v 0 dtsin(q i )), we instead update it (x i + v 0 dtcos(f i ), y i + v 0 dtsin(f i )), where the angle f i is calculated by updating q i with a repulsive displacement f s R/| R| (see [1], for further details and geometrical representation).
Intuitively, this introduces the parameter f s which quantifies the strength of tip-vessel repulsion and represents the tendency of tips to actively move away from regions of high network density (in a typical radius R rep ).However, such repulsive interactions do not change the core dynamics of the first phase of branching morphogenesis.Indeed, the global dynamics of invasive branching morphogenesis by outward propagation of active tips is identical, as is the dynamics of annihilation of tips in the back of the invasive growth front (in 2D, even with repulsion, the geometrical constraints are such that annihilation is always dominant at a given density).The main quantitative di↵erences with repulsion are i) an overall larger network density, as expected from the fact that tips can be partially and locally guided in low density regions prior to termination (see [1] for more details), and ii) an increased local branch alignment driven by repulsion.
Here, we reasoned that this latter e↵ect can be used to constrain the amount of repulsion during lymphatic branching morphogenesis.
Indeed, we computed the nematic order as a proxy for local alignment between individual branch segments in both simulations and experimental data.We quantified the nematic order in network regions of size L ⇥ L and its decay as a function of di↵erent box sizes L. To prevent confounding e↵ects arising from large density fluctuations of the networks, we constrained the analysis for boxes with branch densities comparable to the global density of the network.In each box, we define a local director vector n ⌘ N 1 Â j u j , where u j are individual unit vectors oriented along N discrete branch segments in the box (based on coarse-grained skeletonized coordinates).We then determine the scalar order parameter S ⌘ h2 cos 2 (j j ) 1i averaged over all boxes of size L ⇥ L, where j j is the nematic angle between the unit vectors u j and the local director vector n [5].For P13, we found that the nematic order decayed with a characteristic length estimated by h ' 100 µm from an exponential fit S µ e L/h for small L. Early and late-stage networks P8, P16 and P21 similarly displayed exponential-like decays, with slightly longer decay lengths of h ' 150µm, see Supplementary Figure 4. We reasoned that the faster decay of local nematic alignment at P13 compared with those of P8, P16 and P21 might be indicative of extensive side branching taking place at the onset of the second phase of network growth.We then determined the nematic order in simulations with di↵erent amounts of self-repulsion as quantified by the repulsion strength f s .Interestingly, the decay profile from

Incorporation of branch pruning in the model
To consider the e↵ect of branch pruning (i.e.deletion) in our simulations, we considered that any branch has a finite probability r pr at any given time point to be deleted from the simulations.More specifically, to model the e↵ect of ligand-trap treatment from a specific time point (P11), we ran the simulations with the same parameters as above (Supplementary Figure 4C) and allowed for pruning at di↵erent rates r pr after the second half of the simulations.Only a terminal branch can be pruned, otherwise this would lead to disconnected branch segments.Therefore, for large pruning rates, terminal branch points of lower and lower generation numbers can be selected for pruning, leading to very sparse networks with each subtree only consisting of a a few long branches.Importantly, when calculating the network fluctuation at the end of simulations for di↵erent r pr , we found only marginal e↵ects on the exponent of fluctuations, as expected given that we randomly delete branches from a network already showing giant fluctuations pre-pruning.

Incorporation of side-branching in the model
Given the findings described in the main text, we next sought to incorporate the possibility of side-branching into our model.From a mathematical perspective, this drastically modifies the model, as the state of zero tips (A = 0) is no longer an absorbing state for the dynamics.Intuitively, this means that the local network density is not frozen to a value dictated by the phase of invasive branching morphogenesis, but can continually evolve and increase throughout time.This type of active side-branching is then expected to influence the global topological properties of the growing tree, as explored mathematically in the past [6], and bears similarities with a recent model of branching-delayed random walks that describes salivary gland development [7].
We first consider the simplest case of random side-branching, i.e.where any inactive particle has a random probability r s to become active again (I !A) at any time point.When this occurs, we randomly pick a new growth direction q i perpendicular to the local vessel that the new tip is emanating from.As the process of tip termination still operates, if we have picked a particle that is in a radius R a of other vessels, re-activation of growth/side-branching does not occur.Although adding this new process does not qualitatively change the first phase of tip-driven invasive branching morphogenesis (Figure 2D-G and Supplementary Figure 3G-K), it drastically modifies the second phase by driving a continuous increase of branch number even after the entire ear has been filled up (Supplementary Figure 3I,K).
In particular, if we run this second phase over long times, the branching network converges towards a new maximal density which is dictated by the radius of annihilation R a : all regions of the ear become filled spanning entire regions of the network (with however highly stochastic and broadly distributed sizes, as expected from a process of neutral competition between di↵erent tips during a branching and annihilating random walk).Furthermore, lineage-tracing during the second phase (after the network has reached the edge of the ear) is predicted not to give rise to any large clones, given that all tips are extinct at this point.
In the limit of non-zero side-branching (r s > 0), the picture is slightly di↵erent, as clones can grow even in the second phase when inactive/vessel particles are selected for a side-branching event.However, the picture of wide, bimodal distributions still holds: early lineage tracing is predicted to give rise to a population of very large clones spanning large branching regions (tips being labelled in the first phase of branching) and a population of small clones (either one cell or few cells consisting of local side-branching events).Late lineage tracing on the other hand is predicted to give rise only to this second population of small local clones.
Importantly, this prediction is close to what we observed in our four time points of lineage tracing.Note that we make the typical assumption from the literature that confetti labelled has a time delay of 1.5 day [8,9], and include this when matching the time points between model and data.When lineage tracing is performed during the first phase of branching morphogenesis (P4 and P6), we observe a small subpopulation of broadly distributed clones spanning entire cohesive regions of the network.However, these clones are absent in the later time points (P9 and P12), although we do observe small branch segments being labelled, supporting our side-branching hypothesis (Figure 3A,D).
Although this captures qualitatively the clonal distribution, one limitation with this model is that it gives rise to most clones being of 1 cell size (i.e. cells that were not tip-cells during the simulation).This is however not realistic as proliferation does occur in lymphatic networks (Supplementary Figure 8), as evidenced by the fact that vessel width increases in time (Figure 1A).This increase in vessel width is the reverse of ear growth without lymphatic vessel proliferation, which would be expected to cause vessel thinning.Thus, for the distributions displayed in Figure 3F, we assume a small constant base line proliferation rate of all cells of r d (we choose r d T = 2.5 so that every cell divides a few times during the simulation, in analogy with a division time of a few days).We note that such implementation of division (i.e.locally duplicating particles) does not change the rest of the dynamics, and is computationally costly so we only include it in this part of the analysis.With such division rate, we find that clonal distributions are characterized, at small clone size, by an exponential distribution (as expected from cells with stochastic division patterns [10], as observed in the data, see Figure 3G).
We note that further refinement of the model could be possible.In particular, we assume here that there is a single cell per tip, so that labelling this one cell results in large monoclonal vessels.Although this is broadly consistent with the data, showing indeed that large network regions are monoclonal in P4/P6 inductions, one could consider that "functional tips" are a unit consisting instead of several cells which can move and dynamically exchange position (consistent with experimental observations [11]).As shown previously, this would still give rise to monoclonal network regions after a period of clonal competition, with the length of this period (seen in networks as regions with partially labelled branches prior to full monoclonal conversion) being dependent on the number of functional tip cells [12].

Parameter fitting and model predictions
In this section, we detail how key parameters were fitted in the model based on data, as well as the qualitative and quantitative predictions that the model can subsequently make.

Parameter fitting
We first estimate the dynamics of ear growth, modelled as a half-disk of radius R(t), which grows in time from t = 0 to t = T at the end of simulation (T can also be arbitrarily used to set the time scale of the system).
The initial radius at P4 can be used to set a length scale in the system, so we set it in the simulations to R(0) = 100, and parametrize it according to linear growth (see Figure 1) as R(t) = R 0 (1 + at/T ), where a is the relative expansion at the final time point (P21) compared to the initial time points (P4-P6), which we estimate to a = 1.7.The speed of tip elongation v 0 has a negligible influence in the simulations, as it rescales the dynamics of the system but not its global evolution, so we take v 0 = 1.
On the other hand, the branching probability r b is a key metric as it directly impacts on branch length as well as network density.Experimentally, we observe that average branch length l at P13, when rescaled by ear radius R at that time point, is equal to l/R ⇡ 0.037, which we use in the simulation to fit r b = 0.275 as it gives the right value of l/R at the mid-point of the simulation.As mentioned above, the radius of annihilation has a weak impact on the dynamics, but we estimate it based on the typical distance between terminated tip and the nearest vessel (typically a fraction of the branch length l in our data), leading us to take R a = 3 in our simulation units.For the repulsion radius R rep , we tested di↵erent values for it and found that values much smaller than the typical branch length l result in no changes in the overall network morphology (as repulsion is too local to take place prior to annihilation), while values much larger than branch length were unrealistic given known mechanisms for repulsion (such as local mechanical sensing or local gradient formation typically restricted to a few hundred microns), so we took R rep = 10, i.e. on the order of a branch length.Constraining the value of the repulsion strength f s however is more straightforward, as we found (as described in the section above) that we could use the local nematic order to estimate it, with very large values of f s overestimating the range of branch alignment, while absence of repulsion ( f s ⇡ 0) underestimated it.We thus took f s = 0.2 based on this analysis.
As initial condition, we positioned five tips regularly spaced at the edge of the ear (which itself is delimitated by the curves y = 0, y = Rsinq with q = [0, p]), so that the tip coordinates are (y = 0 and x = 2R 0 /3, R 0 /3, 0, R 0 /3, 2R 0 /3.Each tip points perpendicular to the edge q = p/2, mirroring the experimental conditions of P4 of a few initial trees driving branching morphogenesis (we note that the resulting structure is largely insensitive to this initial condition given the self-organization potential and stochasticity associated with each tip growth).
To directly compare between di↵erent time points in the model and experiment, we simply took P4 as our time point t = 0 and P21 as our time point t = T , and linearly extrapolated in between the two (for instance the middle time point P12.5 corresponding to t = T /2).

Model predictions
With this set of parameters, the model could reproduce a number of important features of the data.Firstly, at the qualitative level, the core assumptions of the model (stochastic branching, tip termination at a small distance from neighbouring vessels) matched well with the non-stereotypical structure of the lymphatic networks in mice, characterised by constant density and front-like propagation of tips during the first phase of branching (as evidenced in P4-P8 in Figure 1).Furthermore, beyond the fact that the model produced networks of constant overall densities after the first phase, the model was characterized by extensive spatial fluctuations (large exponents in the range of 0.6 0.7) in close agreement with the P13 data (Figure 2A,B).At a more quantitative level, we found that i) the local nematic correlation function of vessels was well-described by our simple model of small levels of tip-vessel repulsion during branching morphogenesis (Supplementary Figure 4A,B), ii) the average number of branches per ear after the first phase of branching (P13) was close to the experimentally observed one (Figure 1D, Supplementary Figure 3J) and iii) as described above, clone size distribution showed similar functional shape and time-dependency in model and data (Figure 3F,G).
Perturbation experiments were also found to be in agreement with our theoretical framework.In particular, consistent with the fact that side-branching is the core mechanism responsible for the decrease of giant density fluctuations in the second phase of branching, we found that side-branching inhibition (which could be achieved either by Vegf-c mutant or targeted inhibition of sVEGFR3 from P11 onwards -see Figure 4A-J) results in large fluctuation exponents (a = 0.6 0.7) as predicted in the model with side-branching (first-phase).More specifically, we modelled slightly di↵erently each condition to reflect the timeline of experimental treatment.For Vegf-c mutant, the entire time course of branching morphogenesis is expected to be a↵ected.Therefore, we ran simulations with no side-branching and a branching rate of 25% compared to wild-type (based on Figure 4G), predicting much larger fluctuations as in experiments (Figure 4H-J).For the ligand-trap, targeted inhibition of sVEGR3 was performed at P11, so we keep all parameters as in WT until P11, and then continue to run simulations with no branching (side or tip-branching, leading to regions at the ear edge not yet filled, as seen in data, compare Figure 4D and 4A).Although as shown in Figure 4C, pruning does not significantly a↵ect density fluctuations, we still included it in the model after P11 (to fit the experimental reduction to 40% seen in Figure 4B), as the number of branches already made at P11 in the model would otherwise exceed its experimental value at P21.This also predicted enhanced fluctuations as in the data (Figure 4C, E).
Finally, from a theoretical perspective, we also would predict that increasing branching rate (r b ) via Cpl24 mutants should not qualitatively change the resulting networks, beyond an overall increase in branch density.Indeed, Cpl24 mutant networks were still characterized by local tip termination in the vicinity of vessels, as well as minimal fluctuation exponents (a ⇡ 0.5, consistent with the theoretical argument that this is the minimal value for spatial fluctuation exponent).
3 Local regulation of side-branching: modelling strategy and associated quantification Finally, in this section we provide more information on Figure 5, in particular the di↵erent types of models we considered for the local regulation of side-branching, as well as quantifications done in the data to test these di↵erent models.

Modelling of side-branching initiation
The simplest implementation of side-branching is by taking an equal probability r s for any vessel/inactive particle to become an active tip (if it doesn't sit within distance R a of another vessel).Although such side-branching event can yield to a decrease in giant number fluctuations (as seen in Figure 5A,B), we found quantitatively that this would require a very large increase in the number of branches in the second phase of branching, with typically a two/three-fold increase in branch number required (Figure 5B).
This was in stark contrast to our data, which showed that the reduction in giant density fluctuation exponent a was concomitant with only a 20 30% increase in branch number from P13 to P16/21.We note that this 20 30% increase could be slightly under-estimated due to the possibility that side-branching could have already started between the time point at which the lymphatic network reaches the edge of the ear (P10/11) and the measurement time point of P13.However, this number is still constrained by the fact that networks at P8 already consist of around 700-800 branches despite still low coverage (meaning that the branch number around P10-P11 should be close to P13).
Given previous reports that branching could be induced and/or guided by local factors such as the concentration profile of di↵usible factors or of local hypoxia [13,4,2], we then explored the idea that optimization of space-filling properties could be achieved in a much more parsimonious manner by rendering side-branching events dependent on the local environment of a branch.We thus explored two, non-mutually exclusive, classes of hypotheses on how this could be achieved in vivo: • side-branching rate r s being dependent on the local network density.In this hypothesis, a vessel cell/particle i locally measures the average density r i around it (defined as the number of particles in a given radius R rep divided by the area pR 2 rep ), and has a rate of side-branching inversely related to local density.Many types of numerical implementations are possible for this, including smooth variations of side-branching rate r s as a function of r i , although for the sake of simplicity, we simply take it as binary: below a critical local density r c , side-branches emerge at rate r s while these are prevented (r s = 0) above the critical density threshold.In the numerical implementation of simulations, we thus, as before, iterate over all inactive particles, calculate local density and pick a random number to be compared to r s (r i ), which determines whether branching occurs or not.As in the case of purely constant side-branching, we then pick an orientation q i for the new tip as either perpendicular directions compared to its vessel.
• side-branching rate r s being dependent on the local network density gradient.In this hypothesis, a vessel cell/particle i is not sensitive to absolute local densities, but can read local density gradients.
Density gradients are defined in the exact same way as the implementation of repulsion that we describe above: we calculate the vector R = Â j (r i r j ) which is a weigthed average of all the vectors linking a particle i to its neigbors, and points towards regions of low densities.Vectors with large norms | R| indicate regions where density gradients are large, and we therefore as above make side-branching probability r s (| R|) non-zero only when the repulsion vector is above a threshold | R| c .
In this implementation, the orientation q i of the new tip is taken parallel to R.
Note that in both cases, we assume that the length scale at which cells measure local density (or local density gradients) is R rep , an assumption that could of course be lifted, but is the simplest given the similarity between density-sensing in tip repulsion and density-sensing in side-branching events.As in the case of tip repulsion, taking R rep too large (compared to branch length) would not convey any useful information for side-branching as density would simply converge to its average value, while taking R rep too small would also not give useful information on a cell branching neighborhood.
Interestingly, we found that both implementations of the model drastically increased the speed at which networks converged to minimal fluctuations (characterized by exponents of a ⇡ 0.5).In particular, under both assumptions, low exponents close to 0.5 could be achieved with relatively mild increases in branch number (typically 50 100%.Furthermore, combining both models of density sensing (i.e. for instance by having the magnitude of r s dependent on local density r i , but tip growth direction being picked in the direction of local density gradient R) yielded even better results, with optimal coverage being achieved simply by increasing branch number of 200-300, a number close to the one found in our experimental data (see Figure 1D, 2A,B and 5A,B).

Data analysis and model comparison
To verify these modelling assumptions on side-branching probability being dependent on local density and/or local density gradients, we examined more closely the P13 network dataset, as this is the time point at which we assume from the model that optimization by side-branching would be most important.As discussed in the main text, we could manually locate a number of side-branching events by identifying nascent arrow-headlike branches (Figure 5C).We confined our analysis on the shortest sprouts to prevent potential confounding e↵ects.We performed whole-network reconstructions by skeletonizing the binary images using scikit-image library [14] and the Skan module [15] to extract vector representations of branch segments [1].We then a threshold to identify non-fork-like bifurcation events.
We then located all branching junctions in the skeletonized coordinates of the entire branched networks at P8, P13, P16, and P21, and quantified their vertex configuration factors using a custom-made script.
We determined the angles by examining a small region of ⇠ 8µm surrounding the branching vertex.In addition, we analyzed the vertex configuration factors of manually identified side-branching sprouts at P13.
The fraction of junctions with s 2/3 was highest for the P13 side-branching dataset and lowest for the P8 whole-network data, indicating that symmetric or lateral side-branching configurations are more likely to occur than fork-like tip bifurcation events at later developmental stages (see Supplementary Figure 7D-E).

F)
. B) An image of a P6 mouse head.The white arrow indicates the ventral side of the ear pinna and the blue arrow the tragus that masks the opening to the inner ear.Yellow arrow points at the dorsal side of the ear pinna.C-E) (C) A phase contrast image of a dermis exposed ventral side of a P6 mouse ear pinna.(D) The epifluorescence image of anti-LYVE1 stained ventral ear pinna shows ventral superficial LV networks invading from the base (yellow arrow) and the edge (white arrow).In C and D blue arrow points at the tragus similar to B). E) The ear pinna shown in (C-D) was cleared, imaged with confocal microscope, and depth color-coded.Red LVs grow at the superficial ventral dermis.White arrows point at the connection to the deep dorsal LV network (yellow/green).The original stack is shown in Supplementary Movie 2. Images shown in (C-E) are representative of n=5 ear pinna in 5 mice.Shows a drawing of the lymphatic endothelial sub-trees growing on the superficial ventral dermis in P6 and P8 ear pinna.Each sub-tree is indicated with a unique color.Images are representative of n=5 P6 and n=5 P8 ear pinna and mice.G-I) (G) Shows anti-LYVE1 (red) and anti-VEGFR3 (white) stained P21 ventral ear pinna dermis.Yellow arrows in zoom-in images indicate LYVE1-(pre-)collector segments.Examples of full stacks are shown in Supplementary Movie 7 and 8. H) Graph shows mean proportion +/-SD of the LYVE1+ lymphatic capillary length of all (VEGFR3+) the LVs at P8 (n=5 ear pinna, representing 5 mice), P16 (n=11 ear pinna, representing 9 mice), and P21 (n=10 ear pinna, representing 9 mice).Two-sided LV tree.Yellow arrows point to examples of EdU-positive lymphatic endothelial nuclei.Scale bars for overview and magnified images are 200µm and 50µm, respectively.The shown image is representative of n=5 mouse ear pinna and 5 mice.E) Examples of LYVE1 stained (grey) original image, segmented binary image and the corresponding skeleton used for density fluctuation analysis of the LV-networks at P13 (n=5 ear pinna and mice), P16 (n=8 ear pinna, representing 7 mice) and P21 (n=9 ear pinna representing 8 mice).Scale bars are 500µm.F) Distribution of branch lengths (defined as the distance from two branch points) in reconstructed P13 lymphatic networks, showing broad exponential-like distributions consistent with the model of stochastic branching (black line).Data are shown as mean values +/-SD (n=4 ear pinna in 4 mice).Source data for Supplementary Fig. 2B-C, and F are provided as a Source Data file.average density in the network: left to right in each line for increasing branching rate rb).In all cases, fluctuation exponents are robustly above the minimal value of 0.5 (B, D, F, black line) and consistent with the exponents measured experimentally in wild-type P13 lymphatic networks.G-H) Computational exploration of how different parameters impact the density fluctuation exponent in the presence of side-branching.With intermediate values of side- Fig. 4A-C are provided as a Source Data file.

Supplementary Figure 7 :
, D and G are provided as a Source Data file.Both tip-splitting and side-branching of existing vessels contribute to LV network build-up A-B) LYVE1 staining of a P6 ventral ear pinna dermis.Yellow arrows indicate tip-splitting events and white arrows sprouting of an existing LV.Scale bars are 50µm.A quantification are provided as a Source Data file.A-C) In A) a dot blot shows the mean number +/-SD of A) nascent sprouts of existing LYVE1+ lymphatic capillaries (side branches) (p<0.0001) and B) percent of sprouting lymphatic capillary tips (p=0.0002) in control (n= 4 ear pinna and mice) and ligand trap treated (n= 4 ear pinna and mice) mice, at P13, i.e.only two days after treating mice with VEGF-C/D ligand trap encoding AAV-vector.Samples have been normalized to the size of the analyzed area and then the average of controls has been set to 1. C) LYVE1 staining of P13 ventral ear pinna show invasive tip cell and nascent sprouting of existing LVs (yellow arrows) in control mice, whereas the tips are blunted in VEGF-C ligand trap treated mice (white arrow).D-E) In D) a dot plot shows the mean percentage +/-SD of LYVE1+ lymphatic capillary tips that display a pointy tip cell phenotype upon control (n= 10 ear pinna, representing 6 mice) or VEGF-C ligand trap (n=12 ear pinna representing 6 mice) treatment from P11 to P16 (p<0.0001),i.e. the samples were collected at the time of ongoing LV growth.E) LYVE1 staining of P16 ventral ear pinna show invasive tip cell/sprout phenotype (yellow arrows) in control mice, whereas the tips are blunted in VEGF-C ligand trap treated mice (white arrow).F-H) Dot plot shows mean (+/-SD) F) absolute (p<0.0001) or G) LV area (p<0.0001)normalized (average of controls in each experiment set as 1) number of EdU+ nuclei in the LYVE1+ lymphatic capillaries upon control (n=6 ear pinna, representing 5 mice) or VEGF-C ligand trap (n=11 ear pinna, representing 6 mice) treatment (from P11 to P16).H) A Stacked bar graph shows the proportion of EdU+ LV nuclei in the tip segments and stalk segments of the control (n=4 ear pinna, representing 4 mice) or VEGF-C/D ligand trap treated (n=7 ear pinna, representing 4 mice) mice.I) LYVE1 (red) and collagen IV (white) staining shows an empty basement membrane sleeve (yellow arrow) of an LV in an ear pinna of mice treated with sVEGFR3 from P11 till P21.Images are representative of n=4 ear pinna, representing 3 mice.Scale bars are 20µm.Two-sided Welch's t-test was used for measuring statistical significance in A-B, D, and Supplementary Figure 10: Additional analyses of density-dependent side-branching mechanism at P13 A) Angle between the side-branch initiation vector (see illustration in Fig.5C in the main text) and neighboring branches as a function of sensing radius R for each ear.The color bar represents the relative frequencies of ψ, dotted vertical lines correspond to ψ = ±90 • .B) Ratio of probabilities of neighboring LVs in the "front" (with |ψ| < 45 • ) vs "back" (with |ψ| > 135 • ) of the side branch initiation vector exhibits a conserved length scale of sensing radius R≃200μm for the individual ears.C) Ratio of isotropic densities around side branches ρs to densities ρr around random points on the networks for different values of the sensing radius R. D) Spatial fluctuation exponent in simulations before and after the phase of side-branching (same simulations as Fig. 5A, under the directional and isotropic sensing).Simulation data exhibit the statistical variance of the fluctuation exponent due to the finite system size.Interestingly, the range of the noise attains similar values as in the experimental data, see Fig.2B and S1I.Source data for Supplementary Fig. 10A-D are provided as a Source Data file.